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The next generation of large scale surveys will not only measure cosmological parameters within 
the framework of General Relativity, but will also allow for precision tests of the framework itself. At 
the order of linear perturbations, departures from the growth in the LCDM model can be quantified 
in terms of two functions of time and Fourier number k. We argue that in local theories of gravity, 
in the quasi-static approximation, these functions must be ratios of polynomials in k, with the 
numerator of one function being equal to the denominator of the other. Moreover, the polynomials 
are even and of second degree in practically all viable models considered today. This means that, 
without significant loss of generality, one can use data to constraint only five functions of a single 
variable, instead of two functions of two variables. Furthermore, since the five functions are expected 
to be slowly varying, one can fit them to data in a non-parametric way with the aid of an explicit 
smoothness prior. We discuss practical application of this parametrization to forecasts and fits. 
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I. INTRODUCTION 

Upcoming redshift and weak lensing surveys, such as 
Dark Energy Survey (DES) [1], Euclid [2] and Large Syn- 
optic Survey Telescope (LSST) [3, 4], combined with 
the cosmic microwave background measurements from 
Planck [5] and other cosmological probes, will accurately 
trace the growth of cosmic structures through multiple 
epochs. They will offer the opportunity to test General 
Relativity (GR) by examining the relations between the 
distribution of matter, the gravitational potential and the 
lensing potential on cosmological scales. Such tests may 
yield clues to the physics causing cosmic acceleration or, 
at the very least, extend the range of scales over which 
Einstein's gravity has been validated by experiment. 

To test GR, one can either constrain particular al- 
ternative gravity models, such as the Dvali-Gabadadzc- 
Porrati (DGP) braneworld model [6] or f(R) [7, 8], or 
work within more general parametrized frameworks that 
cover many theories at once and minimize the risk of 
missing potential hints of modified gravity in the data. 
Over the past several years significant effort went into 
developing such frameworks and understanding require- 
ments for their consistency [9-25] . Often, departures from 
the standard cosmological model (LCDM) are quantified 
in terms of arbitrary functions of time and, sometimes, 
scale. These functions cannot be fit to data without fur- 
ther assumptions about their form. 

In this paper we motivate a parametrization that con- 
tains five unknown functions of time only and is general 
enough to cover most viable models of modified grav- 
ity and dark energy proposed so far. Importantly, these 



functions are expected to be slowly varying, hence the ef- 
fective number of degrees of freedom that are fit to data 
can be small. One can avoid assuming a parametric form 
for the five functions and use instead a smoothness prior 
similarly to how it was applied to reconstruction of the 
dark energy equation of state w in [26-28] . 

Observables describing large scale structure are calcu- 
lated using cosmological perturbation theory in Fourier 
space. The relevant variables are the two scalar metric de- 
grees of freedom, e.g., 4> and \l/ in the Newtonian gauge, 
along with the matter density contrast 5 and the mat- 
ter velocity perturbation v. One needs four equations to 
solve for the evolution of these four variables, assuming 
that baryons and dark matter obey the same equations 
at late times. Two equations are provided by the covari- 
ant conservation of matter energy-momentum. The other 
two equations arc supposed to be provided by a theory 
of gravity which prescribes how the metric responds to 
the matter stress-energy. Formally, one can always com- 
plete the system of equations by introducing two func- 
tions p(a, k) and 7(0, k), defined via 1 



k * = -iirfiGa pA , $ = 7* , 



(1) 



where a is the scale factor and A = S + 3aHv/k. They 
are defined in a way that recovers the Poisson and the 
anisotropy equations of LCDM when fi = 7 = 1. There 
are other choices in the literature for the pair of functions 
relating ($, 4<) to (<5, v) that are equivalent to \i and 7. 
A common alternative choice is to use E, defined as 



k 2 ($ + *) = -87rGa 2 E(a, k)pA 



(2) 



in combination with p. As shown in [17, 18, 29, 30], once 
the two functions are given, one has a consistent set of 
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equations that can be incorporated [29-32] into standard 
Boltzmann codes, such as CAMB [33], to calculate the 
observables. In principle, everything that observations 
can tell us about cosmic structure on linear scales can 
be stored as a measurement of fj, and 7 and, if necessary, 
projected onto constraints on specific models. 

But what form should one adopt for these functions 
to fit them to data? In [34-36], a principal component 
analysis (PCA) was performed to forecast the best con- 
strained eigenmodes of /x(a, k) and 7(0, k) for different fu- 
ture surveys, finding that they will measure amplitudes 
of tens (if not hundreds) of them with good accuracy. 
This is encouraging, but it is not clear how many of these 
constrainable eigenmodes are physically interesting. PCA 
alone does not really answer the question of what param- 
eters one should be fitting to data. 

Another concern, which is the main motivation for this 
work, is that an arbitrary relation between two quantities 
in Fourier space, such as those in Eq. (1), does not, in gen- 
eral, imply a local relation between them or their deriva- 
tives in real space. Clearly, the fc-dependence of fi(a, k) 
and 7(0, k) cannot be completely arbitrary if equations 
of motion are obtained from variational principle. 

In this work, we investigate the physically acceptable 
forms of /x(a, k) and 7(0, k) based on considerations of lo- 
cality and general covariance. We show that under rather 
general conditions, and under the quasi-static approxi- 
mation (QSA), they should always have a form of ratios 
of polynomials in k. Furthermore, the numerator of \i is 
set by the denominator of 7. The coefficients inside the 
polynomials are functions of the background quantities 
and can be expected to be slowly varying functions. Tech- 
nically, the number of these time-dependent coefficients 
is infinite if one allows for a completely arbitrary mod- 
ification of GR. However, in models with purely scalar 
extra degrees of freedom, the polynomials are even in 
k and, furthermore, in many viable models considered 
so far in the literature, the polynomials are even and of 
second degree, hence the number of time-dependent co- 
efficients is reduced to five. While this parametrization 
is motivated by the QSA, it allows for departures form 
LCDM on near-horizon scales. 

Baker ct al. [22] have recently investigated the form of 
exact equations of motion for a large variety of modified 
gravity models. They also noted that, under the QSA, 
equations reduce to algebraic relations with even powers 
of k and proposed constraining the time dependence of 
the background dependent coefficients. Instead, we con- 
sider coefficients appearing in fi and 7, which reduces 
the number of free functions in most interesting cases to 
five and makes it easy to use them in existing modified 
Boltzmann codes, such as MGCAMB [29-31]. Amendola 
et al. [24] recently adopted an equivalent five-function 
parametrization to investigate limits of observability of 
modified gravity on linear scales. Their choice was moti- 
vated by results of De Felice et al. [37], who calculated 
/i and 7 in the QSA for the Horndcski [38] class of most 
general second order scalar-tensor theories. We arrive at 



the same form as a particular case of a more general and 
much simpler derivation. 

Working with five arbitrary functions of time may seem 
like a daunting task, but it is much easier than constrain- 
ing two functions of scale and time. Furthermore, the 
functions are known to be slowly varying, which can be 
used as a strong theoretical prior. We outline the prac- 
tical application to forecasts and fits in Sec. IV. Our 
parametrization is useful if one wants to look for depar- 
tures from LCDM without assuming a particular model. 
Clearly, the number of functions can be smaller if one re- 
stricts the range of possibilities. For instance, to describe 
linear perturbations in Brans-Dickc models it is sufficient 
to provide only two functions of the background [39, 40]. 

The remainder of the paper is organized as follows. In 
Sec. II A, we show that, under the QSA, [i and 7 are ra- 
tios of polynomials in k, and the numerator of /x is given 
by the denominator of 7. In II B we point out that for a 
very broad class of viable modified gravity models, the 
polynomials are even and second order in k and, there- 
fore, one needs to specify only five functions of time. In 
Sec. Ill, we examine the assumptions made by the QSA 
and discuss the extent of their applicability. In Sec. IV, 
we outline the procedure for reconstructing \i and 7 from 
data using a smoothness prior applied to the five func- 
tions. We conclude with a discussion in Sec. V. 



II. n AND 7 IN MODIFIED GRAVITY MODELS 
A. The most general case 

Consider a broad class of theories in (3 + 1) dimensions 
with the action defined in terms of a Lagrangian den- 
sity that contains an arbitrary function of geometric in- 
variants R, R a0 R a P, R afSlS R aP7S , AR, R a P\7 a V p R, . . . 
as well as any number of scalar degrees of freedom </>j, 
i = 1, . . . , N (including the longitudinal components of 
vector or tensor degrees of freedom), which can be non- 
minimally coupled to the metric and each other. This 
embraces dark energy models as well as modified gravity 
theories, including effective (3 + 1) dimensional descrip- 
tions of higher dimensional theories. For the moment, 
let us not worry about the existence of ghosts or other 
unphysical properties that such theories may have. At 
this point, we only require invariance of the action un- 
der general coordinate transformations. Let us also make 
an important assumption that there exists a frame in 
which all particles are minimally coupled to the metric, 
so that the matter stress-energy is covariantly conserved. 
For simplicity, we neglect radiation and the differences 
between CDM and baryons. 

Let us now consider the form of equations for linear 
scalar perturbations in the Newtonian gauge, where the 
relevant degrees of freedom of the metric sector are the 
potentials $ and "J defined as 

ds 2 = -(1 + 2^)a 2 dT 2 + (1 - 2$)a 2 dx 2 . (3) 
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Varying the action with respect to the metric tensor gives 
four Einstein equations. The time-time and the time- 
space components can be combined to form the Poisson 
equation, which, to linear order in the perturbations in 
Fourier space, will have the following general form: 



A^ + B$ + C-Sfa = -4irGa 2 pA , 



(4) 



where we assume summation over repeated indices and 
where A, B and C % are linear operators that contain func- 
tions of the background, time derivatives and/or powers 
of k. For instance, for any local theory of gravity, one 
generally has 



A = a-; 



(5) 



where 9™ denotes the time (r) derivative of mth order, 
the highest values of n and m are determined by the 
order of metric derivatives contained in the action, and 
the coefficients a nm are functions of time. Note, that in 
models with only scalar extra degrees of freedom, A will 
contain only even powers of k. This is because odd powers 
can only come from a contraction of spatial derivatives of 
perturbations with spatial indices of the background de- 
pendent coefficients that should vanish for isotropic back- 
grounds such as FRW. Operators B and C l have the same 
form with corresponding coefficients b nm and c l nm . 

Similarly, the traceless space-space Einstein equation 
can generally be written as 



DV + E<S> + F l 5<f> t = 



(G) 



where the zero on the right hand side is due to vanishing 
of the matter anisotropic stress, and the operators D, E 
and F l have the same form as A in Eq. (5). In addition, 
varying the action with respect to each scalar field 4>i, 
and linearizing in the perturbations, will provide equa- 
tions for the corresponding perturbation 5(j>i, that can 
generally be written as 







(7) 



where the operators Hi, Ki and L? also have the form 
given by Eq. (5), with correspondingly renamed coeffi- 
cients. 

Our aim is to find the form of the functions fi(a, k) and 
7(a, k) defined in Eq. (1). Because of the time-derivatives 
in Eqs. (4), (6) and (7), it is impossible to write fi(a 7 k) 
and 7(0, k) in a closed form without solving for the evo- 
lution of the perturbations first. To make progress, let us 
take the quasi static approximation (QSA) in which we 
neglect all time derivatives of ^ and 8 (pi, and delegate 
justifying this approximation to Section III. 

In the QSA, the operators A, B, C\ D, E, F i , H u K l7 
L? become functions, specifically polynomials in k; we 
indicate them with the same letters, removing the hats, 
e.g., we have 



E 



a n0 k n 



(8) 



and similarly for the other functions. As a result, 
Eqs (4), (6) and (7) reduce to a system of linear alge- 
braic equations which we can use to extract fi and 7. 
Defining Ri via 



R^^ 



(9) 



and substituting it, along with $ = 7^, into Eqs. (6) and 
(7), we find 



D + Ej + F l Ri = 
Hi + Ka + LPRj = 



(10) 

(11) 



1 




~E 


F 


-1 


D 


R 




K 


L 




H 



It is convenient to write the solution of these equations 
in the matrix form, 



(12) 



where we introduced a row vector F, a column vector 
K and a square matrix L. We can express the inverse of 
any matrix M as the ratio of its classical adjoint to its 
determinant, M" 1 = adjM/dctM. After some algebra 
we obtain 



dct 



adj 



E F 
K L 



~E 


F 




K 


L 





(E - FL^K) dct L 

E det (L-KE^F) , (13) 

detL -FadjL 
(adjL)A' Ead](L- KE~ l F) ' 

(14) 



Since D, E, F, H , K , L are polynomials, the quantities in 
(13) and (14) are polynomials as well, and, consequently, 
7 and R are fractions of polynomials. Furthermore, in 
their irreducible form, the denominators of 7 and R are 
the same, 



N. 



<?■ R 



Nr 

Q 



where 



A 7 = —D det L 



F(adj L)H , 
N R = L»(adj L)K -E&&i (L - KE~ 1 F)H 
Q =Edet(L- KE^F) . 



(15) 



(16) 
(17) 
(18) 



Also, since the Poisson equation in the QSA has the form 

-4TrGa 2 pA _ k 2 
* p ' 



A + B 1 + C l R t 



(19) 



it follows that fi can be written as an irreducible fraction 
of polynomials with its numerator uniquely determined 
by the denominator of 7 and R, i.e. 



k 2 Q 



AQ + BN 7 + CN E 



(20) 
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As mentioned earlier, in models with pure scalar de- 
grees of freedom, the polynomial functions will contain 
only even powers of k. Furthermore, the k° terms in the 
denominator of \i are negligible in the QSA because the 
magnitude of the corresponding coefficients in the Pois- 
son equation is dependent either on H, time-derivatives 
of <j>i or the first derivative of the scalar field potcntial(s), 
all of which are small in the QSA. This means that the 
k 2 factors in the numerator and denominator of /i will 
cancel. 

Thus, starting from a general covariant action that 
contains an arbitrary function of geometric invariants 
and any number of scalar degrees of freedom (DOF), 
we derived in a model- independent way that, under the 
QSA, the functions 7 and /x are ratios of polynomials in 
k, with the numerator of fi equal to the denominator of 
7- 



B. The subset of viable models 

A parametrization that anticipates a completely arbi- 
trary modification of gravity is impractical, as one cannot 
fit an infinite number of unknown functions to data. Fur- 
thermore, there are good theoretical reasons not to allow 
for arbitrarily high order derivatives or tensorial modes 
in the equations of motion because of the appearance 
of ghost degrees of freedom. Let us therefore consider a 
more representative class of viable modified gravities de- 
scribed by a Lagrangian that contains only one scalar 
DOF obeying second order equations of motion. In this 
case, Eqs (4), (6) and (7) in the QSA reduce to: 

Ak 2 V + Bk 2 ® + Ck 2 6<p = -~4irGa 2 pA , (21) 
£>* + E<£ + F54> = , (22) 
(L + £ifc 2 W = , (23) 



HkH + Kk 2 <& 



where we have made the k dependence explicit, so that A, 
B, . . . , h\ are time dependent coefficients, and we include 
the L coefficient which represents the mass squared of 
the scalar field. 

Following the same steps as in Sec. II A we find in this 
case that [i and 7 are ratios of even polynomials of second 
degree and, as in the general case, the denominator of // 
is the same as the numerator of 7, i.e. 



-DL + (FH - DLi) k 2 



(24) 



' EL + (ELi - KF) k 2 ' 
H = [EL a + (ELi - KF)k 2 ] 1{AE - BD)L Q 

+ [F(BH - AK) + C(KD - EH) 

-1 

+(AE - BD)Li] k 2 1 . (25) 

The above expressions have the same forms as analogous 
expressions derived in [37] for general Horndcski theo- 
ries [38]. Indeed, although we arrived at (24) and (25) 



from general arguments, the subset of viable models to 
which we are restricting coincides with the models in- 
cluded in the Horndcski class, which contains most of 
the viable theories of dark energy and modified gravity. 
The class of theories with a single scalar DOF with a sec- 
ond order equation of motion includes models with ac- 
tions that contain a function f(R, G) of the Ricci scalar 
and the Gauss-Bonnet term, provided the determinant 
of the Hessian is zero, i.e. fimfgg — f%g — 0. Restrict- 
ing to the Lovelock invariants [43] R and Q guarantees 
that no spurious spin-2 ghosts are introduced [42], while 
having a null determinant of the Hessian further en- 
sures that superluminal modes for scalar perturbations 
are avoided [44]. Finally, Horndeski theories include also 
dark energy models such as quintessence and k-essence, 
as well as the covariant Galileon and the 4 dimensional 
effective DGP model in the decoupling limit [41]. 

Without loss of generality, we can rewrite (24) and (25) 
in a more compact way by introducing 5 functions of the 
background {pi(a)}: 



Pi (a) +P2{a)k 2 
l+p 3 {a)k 2 
1+ P3 (a)k 2 



p 4 (o) +p 5 (a)k 2 



(26) 
(27) 



Thus, in the QSA, one can express the perturbed equa- 
tions of motion of a very large class of viable modified 
gravity models in terms of only 5 functions of time. 

Note that even though this ansatz was derived using 
the QSA, it allows for near- and super-horizon modifi- 
cations of gravity: j(a,k — > 0) = pi(a) 7^ 1. Also note 
that, while /z can also deviate from unity on super-horizon 
scales [p, — > p± l {a)), this should not affect any of the ob- 
servables and the super-horizon perturbations will evolve 
consistently with the background expansion [17]. 

Analogous compact forms (26) and (27) were used 
in [24], based on results in [37], to discuss prospects of 
constraining Horndcski theories [38]. We arrive at the 
same forms starting from simpler and more general argu- 
ments that do not require considering the details of the 
Horndeski action. Finally, the form of 7 and [i in (26) 
and (27) resembles the parametrization introduced by 
Bertschinger and Zukin (BZ) in [12], however, there are 
some important differences. In the BZ parametrization, /x 
and 7 tend to 1 for k — >• 0, effectively reducing the theory 
to GR in that limit. Furthermore, BZ does not set the 
numerator of \i equal to the denominator of 7, which we 
find instead to be a general feature. Finally, they fix the 
time dependence of the coefficients in the k 2 expansion 
to a power law, while we leave it general. From a theo- 
retical point of view, not all time dependencies can be 
described as power laws; however, the differences might 
be undetectable depending, in part, on the range of red- 
shifts probed by the experiment. 
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III. THE QUASI-STATIC APPROXIMATION 

As mentioned earlier, closed form expressions for 
/J,(a,k) and 7(0, k) may not exist in a general gravity 
theory unless one adopts the QSA, in which the relations 
between the metric potentials and the matter perturba- 
tion become algebraic. Without the QSA, one needs to 
solve the differential equations in order to define \i and 
7, making them dependent on the initial conditions. 

There are two different assumptions involved in what 
is commonly known as the QSA: (f ) the relative small- 
ness of the time derivatives of metric perturbations com- 
pared to their space derivatives, and (2) the sub-horizon 
approximation, k/aH ^> 1. In LCDM, the second as- 
sumption automatically implies the first — the perturbed 
quantities evolve on time scales comparable to the ex- 
pansion rate, thus their time derivatives become com- 
parable to space derivatives only for perturbations on 
near-Hubble scales. In alternative gravity models with 
additional degrees of freedom, the two assumptions need 
not imply each other, so let us separately consider their 
effects on the applicability of our paramctrization. 



A. Neglecting time-derivatives 

In scalar-tensor models of modified gravity, there can 
be rapid oscillations of the metric potentials on top 
of their slow evolution which can make their time- 
derivatives large. However, the terms that set the am- 
plitude and the frequency of such oscillations are propor- 
tional to the strength of the coupling and the range of the 
extra scalar degree of freedom. Because of the generally 
tight constraints on fifth forces, one typically finds that 
the oscillations are undetectable. 

To illustrate the point, let us consider models in which 
the Einstein-Hilbert part of the action is given by R + 
f(R), and the Poisson equation has the form [45] 



fc 2 * + fc 2<% 



3 

2F + 2 



F ■ 



(28) 



where /a = df /dR = F — 1 is the "scalaron" degree of 
freedom, T-L = a~ 1 da/dr, and the non-quasi-static terms 
are collected inside the square brackets. On sub-horizon 
scales (when k <C aH), the first term inside the square 
brackets is negligible compared to the k 2 5fn term. But 
on large scales it is still smaller than other terms because 
d.fn = JrrSR, and fun must be small. The latter is re- 
quired for the Chameleon mechanism [46] to screen the 
fifth force inside our solar system — the values of 
and fun must be small [47]. The second term inside the 
square brackets is large only if <f> is large. However, the 
evolution of $ in the heavy scalaron (small Jrr) limit 
is practically the same as in GR, except for additional 
oscillations with a tiny amplitude set by /&r [48]. Even 



if the oscillations had a larger amplitude, they would be 
difficult to detect because of their high frequency set by 
Irr- Furthermore, there are no oscillations in the lensing 
potential $ + ty, hence there can be no signal in the In- 
tegrated Sachs- Wolfe effect that constrains the near- and 
super-horizon evolution of perturbations. 

We are not aware of a theory in which oscillations in 
extra DOFs are observable for the range of parameters 
that has not already been ruled out. Thus, it is reason- 
able to adopt a theoretical prior that ignores rapid time- 
variations of gravitational degrees of freedom until we 
find an example of a viable theory in which they are ob- 
servable. Finding such an example may warrant an ap- 
propriate extension of Eqs. (26) and (27). 



B. The sub- horizon approximation 

By ignoring the time derivatives in the modified equa- 
tions we are neglecting not just the rapid oscillations in 
metric perturbations but also the slowly varying signa- 
tures of modified gravity. This, in the absence of addi- 
tional information about the model, can only be justified 
in the k/aH 3> 1 limit. 

Before addressing the significance of near- horizon mod- 
ifications of gravity, let us make an important point: the 
implementation of our paramctrization in the equations 
of motion does not assume the QSA. In the LCDM limit, 
when fi — 7 — 1, we recover the exact equations of 
GR, while the parametrization allows for departures from 
(i = 7 = 1 on all scales. Also, we do not suggest that one 
should ignore the relativistic effects when calculating the 
observables [49-55] . The implementation of near- horizon 
and other relativistic effects [51] in Boltzmann codes like 
CAMB is unaffected by the use of the (7, fj,) parametriza- 
tion — only the Einstein equations are modified, while 
Boltzmann equations and the expressions for the observ- 
able quantities remain the same as before. 

The two relevant questions about the validity of taking 
the QSA limit in deriving the form of our parametrization 
are: (1) How observable departures from the LCDM pre- 
diction can be on near horizon scales in viable modified 
gravity models? and (2) To what extent our paramctriza- 
tion would bias such a potentially observable signature? 

As far as we know, there is no example of a the- 
oretically motivated model that is not ruled out and 
in which departures from LCDM on near-horizon scales 
have been shown to be detectable [59], although there is 
clearly more room for investigation. Part of the reason is 
that cosmic variance limits the statistical significance of 
any inference on large scales. Multi-tracer technique pro- 
posed in [56-58] can remove this limitation to some ex- 
tent. However, in order for the modified gravity signal on 
large scales to be detectable, it has to be sufficiently pro- 
nounced while still keeping the model in agreement with 
other constraints. While one can design specific models 
providing such an example [59], they cannot be consid- 
ered representative. 
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The second question - - the extent to which our 
parametrization would bias a potentially important sig- 
nature on near-horizon scales — can only be answered by 
considering particular solutions of specific models. But 
first one has to find a viable model in which such signa- 
tures are observable at all. 

A conservative way to use the parametrization in 
Eqs. (26) and (27) would be to separately fit to a subset 
of data corresponding to clustering on sub-horizon scales. 
Then, if a departure from LCDM is seen, one would have 
a clear idea about which scales contribute the most to 
the signal, and whether it is appropriate to interpret it 
under the QSA. 



IV. PRACTICAL APPLICATION TO 
FORECASTS, CONSTRAINTS AND 
RECONSTRUCTIONS 

It is impossible to constrain a function without assum- 
ing something about its form. One possibility is to pick 
a particular functional form, such as a power law depen- 
dence, similar to how it is done in the BZ parametriza- 
tion [12]. Since the five functions in Eqs. (26) and (27) 
are expected to be slowly varying, this need not be a bad 
approximation if the data only probes a limited range 
of redshifts. But it is unlikely that a single power law 
will capture the evolution over a wide range of epochs, 
which is what we can expect from surveys like Euclid, 
SKA or LSST in combination with CMB and other data. 
Let us instead explore a way to constrain /j, and 7 in 
a non-parametric way that still takes into account their 
smoothness. 

Recent Refs. [27, 28] proposed a transparent Bayesian 
framework for constraining the dark energy equation of 
state w(a) based on adopting an explicit smoothness 
prior. The prior is defined via a correlation function that 
correlates values of w at neighbouring points in a. This 
framework can be applied to any unknown function (or 
functions) expected to be smooth from theoretical consid- 
erations. Let us outline how it can be applied to {pi(a)} 
in Eqs. (26) and (27) for the purpose of forecasting future 
constraints on 7 and /1, as well as for fitting to real data. 



A. Application to forecasting 

As a starting point, one can discretize the functions 
into finite numbers of bins. The binning can be imple- 
mented using a smooth function, such as a hyperbolic 
tangent, to avoid infinite derivatives at the edges, and 
the number of bins can always be taken to be suffi- 
ciently large to achieve convergence 2 . The binned values 



2 One of the advantages of the smoothness prior approach of [26- 
28] is that it eliminates the dependence on binning. 



of {pi(a)} can be substituted into Eqs. (26) and (27) to 
find fj, and 7 which are used as input in a modified Boltz- 
mann code such as MGCAMB [31]. Along with providing 
other cosmological parameters, this is sufficient for cal- 
culating all types of cosmological observables. 

In the simplest approach to forecasting, one assumes a 
Gaussian shape of the parameter likelihood surface with 
a peak corresponding to a fiducial model, and proceeds 
to calculate the Fisher matrix from derivatives of ob- 
servables with respect to model parameters. One then 
inverts it to obtain an estimate of the total covariancc 
matrix. Because of the large number of highly correlated 
parameters, considering a constraint on any single bin is 
meaningless. One can instead use the principal compo- 
nent analysis (PCA) [60] to see which independent lin- 
ear combinations of bins will be best constrained by a 
given experiment. This is accomplished by diagonalizing 
the corner of the covariance matrix corresponding to the 
bins of the five functions. 

What can one do with the information obtained from 
the PCA forecast for {pi(a)}7 There will be a strong 
degeneracy between the five functions, which means one 
should look at independent linear combinations of bins 
of all five functions. The number of such well-constrained 
combined eigenmodes should give us a measure of how 
many physically relevant independent degrees of freedom 
one can measure about /i and 7. 

In the absence of a theoretical prior, all eigenmodes of 
{pi(a)} carry some information. How should one decide 
which modes are informative and which are not? This is 
where one can use the knowledge about the slowly vary- 
ing nature of the functions and introduce a smoothness 
prior [26-28] . The prior comes in a form of a non-diagonal 
matrix C pmr that one adds to the original covariance ma- 
trix and which introduces additional correlations between 
bins of the five functions. One then considers an eigen- 
mode to be informative if it is unaffected by the prior, 
i.e. if it is the same before and after addition of the prior 
covariance. 

To be more explicit, let us assume that each function 
Pi{a) is binned into N bins in the scale factor, a\, . . . , ajv, 
and label the bins so that they form a single vector p a 
with a = 1, . • • , 5N. For example, one can define 

pi =p 1 (a 1 ), p 2 =Pi(a 2 ),...,PN =Pi(cin), 

Pn+i =£2(01 ),.-., P5N =Pb (cln)- (29) 

Let C^ ta be the corner of the covariance matrix corre- 
sponding to {p a } that we have obtained earlier by invert- 
ing the total Fisher matrix. According to Bayes' theorem, 
the posterior probability distribution for the parameters 
{p a } is a product of the likelihood and the prior prob- 
ability. For Gaussian probabilities, this implies that the 
net covariance matrix (i.e. corresponding to the posterior 
probability) is a sum of C* data and the prior covariance 
C prlor . To construct the latter, one can follow the pre- 
scription in [26, 27] and start with specifying a correla- 
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tion function 

(foCo) -pf^ipiia') -p^ m )) = - a'\) (30) 

where p l ^ dm are the constant values of Pi(a) in LCDM 
(pi = P4 = 1 and p2 = P3 — P5 = 0), and where the form 
of the functions $,^\\a — a'\) is chosen so that 

£ w (Aa) £ w (0) when Aa = |a - a'| < a c 
£ (i) (Aa) ->■ when Aa » a c , 

where a c is the correlation length and £^(0) is a posi- 
tive constant 3 . One has to specify the functional form of 
^ (Aa) and we refer the reader to [27] for an extended 
discussion of different choices. The choice does not make 
a big difference in practice. Using ^''(Ao), one can cal- 
culate the prior covariance matrix for the binned Pi(dj) 
via 

where da is the width of a bin in a. One can define such an 
N x N matrix for each of the five functions and combine 
them to form a block diagonal 5N x 5N matrix C prlor for 
the parameters p a such that 

C prior = diag[C (1) ,...,C (5) ] . (32) 

This prior assumes that the five functions are indepen- 
dent of each other, which is true in the most general case 
but not in many specific models. For example, in f(R) 
only one function is independent, while in more general 
Brans-Dicke models there are two. Such additional re- 
strictions can be implemented, if desired, by adjusting 
the form of C" prior . 

Having constructed the prior matrix, one then com- 
pares the eigenmodes of C data to the eigenmodes of 
£(data_|_£<prioi^ rp ng e jg enm0( j es that are common to both, 

i.e. those that survive, can be considered to be informa- 
tive with respect to that prior. Due to the nature of the 
prior, one expects that slowly varying eigenmodes that 
are best constrained by data will have a higher chance 
to survive, while the high frequency modes will be sup- 
pressed. Naturally, the outcome of this comparison de- 
pends on the parameters of the prior, a c and £^'(0). 
Their choice should, in principle, come from our theo- 
retical prejudice. In practice, they can be tuned so that 
eigenmodes with variations on time scales comparable to 
Hubble time (or another time-scale that is theoretically 
motivated) survive. Thus, a PCA forecast is a key step 
for tuning the prior that can later be used in fitting to 
data, as we discuss next. 



In certain cases, it may be more appropriate to specify the cor- 
relation between points in log a rather than in a. 



B. Fitting to real data 

In a Fisher forecast, there is no limitation on making 
the number of bins N as large as needed to achieve a 
convergence of well-constrained eigenmodes to the con- 
tinuous limit. But this is not the case when fitting to 
real data, which involves searching for the maximum of 
a multi-parameter likelihood surface. The search stalls if 
the parameter space contains flat directions correspond- 
ing to nearly degenerate combinations of parameters, 
which is guaranteed to be the case when the number of 
bins is large. Fitting only the few best constrained eigen- 
modes amounts to a rather strong assumption that the 
amplitudes of the poorly constrained modes are known to 
be exactly zero, which amounts to adopting a strong, yet 
somewhat obscure, prior. Instead, in [27, 28] it was sug- 
gested to use the explicit smoothness prior described in 
the preceding subsection to aid the convergence of Monte 
Carlo Markov chains (MCMC). This is achieved in prac- 
tice by adding a term 

X*nor = (P - p lcdm ) T [CP riOT ]- 1 (p - p lcdm ) (33) 

to the Xdata m MCMC. The number of bins that one fits 
to data need not be very large. One typically needs a 
couple of bins per effective correlation length set by the 
prior. It remains to be shown for specific experiments 
how strong the prior needs to be in order for MCMC to 
converge. Once MCMC has converged, one can quantify 
the statistical significance of the detection of a departure 
from LCDM from the improvement in the Xdata- *-* ne 
can also compute the evidence for the best fit model and 
the Bayes' factors, since the prior probability is explicitly 
known. An explicit illustration of such a calculation for 
the case of w(a) can be found in [28]. 

We do not expect the reconstructed shapes of the 
individual functions Pi(a) to be highly informative be- 
cause observables will only constrain their combinations 
and degeneracies between parameters {p a } will make 
marginalized errors on them large. One could instead use 
Eqs. (26) and (27) to visualize reconstructions of \i and 
7, or fi and S, as surfaces in the (a, k) space. 



V. SUMMARY 

In this paper we have motivated a parametric form for 
the modified growth functions \i and 7 that fixes their 
scale dependence to a ratio of polynomials in k and have 
shown that generally the denominator of 7 is equal to 
the numerator of fi. We arrive at this form by taking 
the quasi-static approximation (QSA) in the equations 
for scalar perturbations derived from a covariant action 
that allows for modifications of gravity and any number 
of scalar degrees of freedom. We examine the impact of 
assuming the QSA in our derivation and conclude that, 
until a viable counterexample is found, the non-quasi- 
static effects of modified gravity can be assumed to be 
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negligible. We nevertheless note that the final form of 
our parametrization allows for a detection of some non- 
quasi-static signatures, but not necessarily for the most 
general ones in that regime. 

We further argue that for most of the viable modifi- 
cations of gravity discussed in the literature, the polyno- 
mials in k are even and of second degree, effectively re- 
ducing the number of time-dependent coefficients to five. 
Since these coefficients are functions of the background 
variables only, they can be safely assumed to be slowly 
varying. This justifies using an explicit smoothness prior 
on their shape when fitting them to data, similarly to the 
Bayesian framework developed for the reconstruction of 



w(a) in [27, 28] . 

A forecast of future reconstructions of fi and 7 based 
on this approach is currently in progress [61]. 
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